Supplementary 3: Modelling TOC concentration and forecast in coastal drainage basins
Camille M. Crapart
24 06 2022
library(dplyr)
library(tibble)
library(DBI)
library(sf)
library(gimms)
library(ncdf4)
library(rgdal)
library(stringr)
library(ggplot2)
library(RColorBrewer)
library(colorspace)
library(cowplot)
library(png)
library(spatialreg)
library(kableExtra)
library(reshape2)
library(grid)
library(gridExtra)
library(ggpubr)
source("f_plotspatial.R")
memory.limit(size = 160000)
## [1] Inf
library(captioner)
fig_nums <- captioner()
table_nums <- captioner(prefix = "Table")1 Coastal drainage basins data
1.1 Connections to database
The coastal drainage basins polygons were computed with a seamless digital elevation model for Fennoscandia (Finstad 2017). The model resulted in more than 200 000 coastal drainage basins, most of them being extremely small polygons covering the coastal regions. In order to simplify the model and reduce the amount of data, we keep only the 0,001 % of the coastal drainage basins that are the biggest. They account to more than 93% of the total surface covered by the initial amount of polygons.
con <- dbConnect(RPostgreSQL::PostgreSQL(),user = "camille.crapart", password = "camille",host = "vm-srv-wallace.vm.ntnu.no", dbname = "nofa")
query <- paste("SELECT gid,geom FROM ","\"Hydrography\"",".waterregions_dem_10m_nosefi",sep = "")
query.2 <- paste("SELECT * FROM ","\"Hydrography\"",".Fenoscandia_LakeRegister",sep = "")
water.regions <- st_read(con,query = query)
lake.register <- st_read(con,query = query.2)
saveRDS(water.regions,"water.regions.Rdata")water.regions <- readRDS("water.regions.Rdata")
ggplot()+geom_sf(data=water.regions,aes(fill = gid))water.regions <- readRDS("water.regions.Rdata")
water.regions$area <- st_area(water.regions) %>% as.numeric()
area1 <- sum(water.regions$area)
wr.quantile <- quantile(water.regions$area, probs = c(0.99,0.995,0.999))
wr.sf <- water.regions %>% filter(area > wr.quantile[3])
area2 <- sum(wr.sf$area)
ratio <- area2/area1*100
map.wr <- ggplot(wr.sf) + geom_sf(aes(fill=gid,col=gid))
ggsave(plot = map.wr, filename = "WR/wr.sf.png",dpi = "retina", width = 10, height = 15, units = "cm")
saveRDS(wr.sf,"WR/wr.sf.rds")
st_write(wr.sf, "wr.sf.shp")1.2 NDVI
NDVI values are extracted from the GIMMS NDVI3g dataset (The National Center for Atmospheric Research 2018), stored on http://poles.tpdc.ac.cn/en/data/9775f2b4-7370-4e5e-a537-3482c9a83d88/, in the “ecocast” dataset and accessed via the “gimms” package (Detsch 2021) on Rstudio.
Data is taken bi-monthly (Tucker et al. 2005). Raster layers of all slices are downloaded using the gimmsdownload function, then monthly composite are calculated using the monthlyComposites function. The maximum NDVI value is taken for each pixel. Afterwards the NDVI values for each catchment polygon are extracted using the extract function from the raster package (Hijmans et al. 2018).
The values stored in the ecocast dataset are composed of the NDVI value and a flag value indicating the goodness of the data (Detsch 2021). All the NDVI values extracted, corresponding to the values for the studied catchments, had a flag of 1, indicating a good value. The NDVI value was retrieved from the stored value using the formula \(NDVI = floor(ndvi3g/10)/1000\) (Detsch 2021), and then the mean of the 3 summer values (June, July and August) was computed for each catchment.
ndvi.max <- readRDS("ndvi.max.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
summer.mean <- raster::stack(ndvi.max[[6]],ndvi.max[[7]],ndvi.max[[8]]) %>% mean()
summer.ndvi <- reclassify(summer.mean,cbind(-Inf, 0, NA), right=FALSE)
summer.scandinavia <- raster::crop(summer.ndvi,c(0,35,55,73))
writeRaster(summer.mean,"ndvi1994.tif",overwrite = T)
ndvi.wr <- raster::extract(summer.ndvi,wr.sf,fun = mean, df = T, sp = T, na.rm = T)
ndvi.wr$area <- NA
ndvi.wr$ndvi <- (floor(ndvi.wr$layer)/10)/1000
saveRDS(ndvi.wr, "WR/ndvi.wr.rds")
knitr::include_graphics("WR/wr.map.NDVI.png")
Figure 1: Summer NDVI in 1994 in coastal drainage basins
1.3 Runoff
The runoff data was downloaded from the CORDEX database: https://portal.enes.org/data/enes-model-data/cordex. The data is available at https://esg-dn1.nsc.liu.se/projects/esgf-liu/.
The names of the variables in the CORDEX project follow the CF Metadata convention: http://cfconventions.org/
THe mean of the 30 years between 1970 and 2000 is used (p11 on EURO-CORDEX guidelines https://www.euro-cordex.net/imperia/md/content/csc/cordex/guidance_for_euro-cordex_climate_projections_data_use__2021-02_1_.pdf). This time period matches the temperature and precipitation data provided by WorldClim.
runoff.raster <- "CORDEX/mrros_Lmon_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_185001-201412.nc"
runoff.nc <- nc_open(runoff.raster)
runoff.stack <- raster::stack(runoff.raster, bands = c(1441:1812))
runoff.mean <- runoff.stack %>% mean()
wr.sf <- readRDS("WR/wr.sf.rds")
runoff.wr95 <- raster::extract(runoff.mean, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(runoff.wr95) <- c("gid","area","runoff")
saveRDS(runoff.wr95,"WR/runoff.wr95.rds")knitr::include_graphics("WR/wr.map.Runoff.png")
Figure 2: Runoff in 1995 in coastal drainage basins (mm/y)
1.4 Land Cover
The Land Cover data was downloaded from https://land.copernicus.eu/pan-european/corine-land-cover/lcc-2000-2006. The previous version (1995-2000) excludes Norway so the 2000 version was preferred.
The land cover data was downloaded as a raster (.tiff file), which included the legend recording the code for each land use. The raster and legend were assembled in QGIS before being imported as a shapefile in R. The codes (corresponding to the line in the extracted dataframe) of the categories of interest were:
Arable land:
- 12: non-irrigated arable land
- 16: fruit trees and berries plantations
- 18: pastures
- 19: annual crops associated with permanent crops
- 20: Complex cultivation patterns
Bogs:
- 36: peat bogs
- 35: inland marshes, excluded because all zeros in the studied area.
Forest:
- 23: Broad-leaved forest
- 24: Coniferous forest
- 25: Mixed forest
Bare:
- 31: Bare rocks
- 32: Sparsely vegetated areas
- 33: Burnt areas
wr.sf <- readRDS("WR/wr.sf.rds")
corine <- raster::raster("Corine_Land_Cover/U2006_CLC2000_V2020_20u1.tif")
wr.corine <- st_transform(wr.sf, crs(corine))
extract_corine <- function(x){
clc <- extract(corine,x)
saveRDS(clc,file = paste("WR/corine/clc",x$gid,"rds",sep="."))
}
wr.list <- split(wr.corine,seq(nrow(wr.corine)))
lapply(wr.list[1143:1440],extract_corine)
#clc.test <- readRDS("WR/clc/clc.16711649.rds")
clc.files <- list.files(path = "WR/corine/", pattern = "clc.*.rds", full.names = T)
list.clc <- sapply(clc.files,readRDS)
list.clc.wr <- sapply(list.clc, as.integer)
names.list <- str_extract(clc.files,"\\d+")
names(list.clc.wr) <- names.list
saveRDS(list.clc.wr,"WR/corine/list.clc.wr.rds")list.clc.wr <- readRDS("WR/corine/list.clc.wr.rds")
clc.tab.wr <- sapply(list.clc.wr, function(x) tabulate(x,45)) %>% t()
clc.tab.area.wr <- clc.tab.wr*prod(res(corine))
catchment.area <- rowSums(clc.tab.area.wr)
clc.tab.prop.wr <- sweep(clc.tab.area.wr,1,catchment.area, FUN = "/")
saveRDS(clc.tab.prop.wr,"WR/corine/clc.tab.prop.wr.rds")clc.tab.prop.wr <- readRDS("WR/corine/clc.tab.prop.wr.rds") %>% as.data.frame() %>% tibble::rownames_to_column(var = "gid")
bogs <- clc.tab.prop.wr %>% dplyr::select(c("gid","V36")) %>% setNames(c("gid","bogs"))
saveRDS(bogs,"WR/corine/bogs.wr.rds")
arable <- data.frame( gid = clc.tab.prop.wr$gid, arable = rowSums(select(clc.tab.prop.wr, c("V12","V16","V18","V19","V20"))))
saveRDS(arable,"WR/corine/arable.wr.rds")
forest <- data.frame( gid = clc.tab.prop.wr$gid, forest = rowSums(select(clc.tab.prop.wr, c("V23","V24","V25"))))
saveRDS(forest,"WR/corine/forest.wr.rds")
bare <- data.frame( gid = clc.tab.prop.wr$gid, bare = rowSums(select(clc.tab.prop.wr, c("V31","V32","V33"))))
saveRDS(bare,"WR/corine/bare.wr.rds")
water <- data.frame( gid = clc.tab.prop.wr$gid, water = rowSums(select(clc.tab.prop.wr, c("V40","V41","V42","V43","V44"))))
saveRDS(water,"WR/corine/water.wr.rds")list.grob <- c("WR/wr.map.Bog.png", "WR/wr.map.Arable.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow = 1, ncol = 2, labels = c("a)","b)"), label_size = 10)
Figure 3: Proportion of a) bogs and b) arable land in 2000 in coastal drainage basins
1.5 Temperature and precipitation
Temperature and precipitation are downloaded from the WorldClim database: https://worldclim.org/data/index.html
Temperature in historical Worldclim is in °C * 10 (to reduce file size). Precipitation are average monthly precipitation in mm. Here we use the mean monthly precipitation and temperature for the period 1970-2000.
wr.sf <- readRDS("WR/wr.sf.rds")
bio.fennoscandia <- raster::getData(name = "worldclim", var = "bio", res = 2.5)
temp.wr <- raster::extract(bio.fennoscandia$bio1, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(temp.wr) <- c("gid","area","temp")
temp.wr$temp <- temp.wr$temp / 10
saveRDS(temp.wr,"WR/temp.wr.rds")
#temp.wr <- readRDS("WR/temp.wr.rds")
prec.wr <- raster::extract(bio.fennoscandia$bio12, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(prec.wr) <- c("gid","area","prec")
saveRDS(prec.wr,"WR/prec.wr.rds")
#prec.wr <- readRDS("WR/prec.wr.rds")
tp <- cbind(wr.sf[,1],temp.wr[,3]) %>% cbind(prec.wr[,3])
saveRDS(tp,"WR/tp.rds")list.grob <- c("WR/wr.map.Temp.png", "WR/wr.map.Precip.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow= 1, ncol = 2, labels = c("a)", "b)"), label_size = 10)
Figure 4: Map of a) average temperature (°C) and b) average precipitation (mm) in 1995 in coastal drainage basins
1.6 Nitrogen deposition (TNdep)
The nitrogen deposition data was extracted from the EMEP database (https://emep.int/mscw/mscw_moddata.html#Comp). The data from 2000 was used as it is the first model available.
N-deposition (NOx and NH3) is the sum of:
- Dry deposition of oxidized nitrogen per m2 grid DDEP_OXN_m2Grid, in mg/m2
- Wet deposition of oxidized nitrogen WDEP_OXN, in mg/m2
- Dry deposition of oxidized nitrogen per m2 grid DDEP_RDN_m2Grid, in mg/m2
- Wet deposition of reduced nitrogen WDEP_RDN, in mg/m2
library(ncdf4)
EMEP_file <- "EMEP/EMEP01_rv4.42_year.2000met_2000emis_rep2021.nc"
wr.sf <- readRDS("WR/wr.sf.rds")
woxn <- raster(EMEP_file, varname = "WDEP_OXN")
doxn <- raster(EMEP_file, varname = "DDEP_OXN_m2Grid")
wrdn <- raster(EMEP_file, varname = "WDEP_RDN")
drdn <- raster(EMEP_file, varname = "DDEP_RDN_m2Grid")
woxn_wr <- extract(woxn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(woxn_wr) <- c("gid","area","woxn")
saveRDS(woxn_wr,"woxn_wr.rds")
doxn_wr <- extract(doxn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(doxn_wr) <- c("gid","area","doxn")
saveRDS(doxn_wr,"doxn_wr.rds")
wrdn_wr <- extract(wrdn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(wrdn_wr) <- c("gid","area","wrdn")
saveRDS(wrdn_wr,"wrdn_wr.rds")
drdn_wr <- extract(drdn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(drdn_wr) <- c("gid","area","drdn")
saveRDS(drdn_wr,"drdn_wr.rds")
ndep.wr <- merge(woxn_wr,doxn_wr, by = "gid") %>% merge(wrdn_wr, by = "gid") %>% merge(drdn_wr, by = "gid")
ndep.wr$tndep <- rowSums(ndep.wr@data[,c(3:6)])
saveRDS(ndep.wr,"WR/ndep.wr.rds")knitr::include_graphics("WR/wr.map.TNdep.png")
Figure 5: Total nitrogen deposition in 1995 in coastal drainage basins (mg/m2)
1.7 Elevation
The altitude data was accessed through the “getData” function, and comes from the NASA digital elevation model (SRTM 90 m). The slope was then computed using the “terrain” function from the raster package.
wr.sf <- readRDS("WR/wr.sf.rds")
alt.nor <- raster::getData("alt",country = "NOR", mask = T)
alt.swe <- raster::getData("alt",country = "SWE", mask = T)
alt.fin <- raster::getData("alt",country = "FIN", mask = T)
alt.list <- list(alt.nor,alt.swe,alt.fin)
alt <- do.call(raster::merge,alt.list)
alt.wr <- raster::extract(alt, dplyr::select(wr.sf, c("gid","geom")), sp = T, fun = mean, df = T, na.rm = T)
names(alt.wr) <- c("gid","alt")
saveRDS(alt.wr,"WR/alt.wr.rds")
slope.nor <- terrain(alt.nor, unit = "degrees")
slope.swe <- terrain(alt.swe, unit = "degrees")
slope.fin <- terrain(alt.fin, unit = "degrees")
slope.list <- list(slope.nor,slope.swe,slope.fin)
slope <- do.call(raster::merge,slope.list)
slope.nona <- reclassify(slope,cbind(NA,100), right = NA)
slope.wr <- raster::extract(slope,wr.sf, sp = T, df = T, fun = mean, na.rm = T)
names(slope.wr) <- c("gid","area","slope")
saveRDS(slope.wr, "WR/slope.wr.rds")1.8 Gathers dataset of coastal drainage basins
The data was gathered and merged into a single dataframe, based on the catchment identifyer (“gid”).
ndvi.wr <- readRDS("WR/ndvi.wr.rds")
ndvi.wr$area <- NULL
runoff.wr95 <- readRDS("WR/runoff.wr95.rds")
runoff.wr95$area <- NULL
bogs.wr <- readRDS("WR/corine/bogs.wr.rds")
forest.wr <- readRDS("WR/corine/forest.wr.rds")
arable.wr <- readRDS("WR/corine/arable.wr.rds")
bare.wr <- readRDS("WR/corine/bare.wr.rds")
water.wr <- readRDS("WR/corine/water.wr.rds")
#clc.wr.sum <- readRDS("WR/clc/clc.wr.sum.Rdata")
alt.wr <- readRDS("WR/alt.wr.rds")
slope.wr <- readRDS("WR/slope.wr.rds")
slope.wr$area <- NULL
tp <- readRDS("WR/tp.rds")
ndep.wr <- readRDS("WR/ndep.wr.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
all.wr.sf <- merge(wr.sf,runoff.wr95,by = "gid") %>% merge(ndvi.wr,by="gid") %>% merge(bogs.wr,by = "gid") %>% merge(forest.wr, by = "gid") %>% merge(arable.wr, by = "gid") %>% merge(bare.wr, by = "gid") %>% merge(water.wr, by = "gid") %>% merge(st_drop_geometry(tp),by = "gid") %>% merge(alt.wr, by = "gid") %>% merge(slope.wr, by = "gid") %>% merge(ndep.wr, by = "gid")
saveRDS(all.wr.sf, "WR/all.wr.sf.rds")The data was then cleaned and the centroid of each coastal drainage basin was computed.
all.wr.sf <- readRDS("WR/all.wr.sf.rds")
wr.sf.95 <- all.wr.sf %>% dplyr::filter(ndvi != 0) %>% filter(runoff > 0) %>% filter(is.na(bogs) == F) %>% filter(is.na(alt) == F) %>% filter(is.na(slope) == F)
wr.sf.95 <- wr.sf.95[!duplicated(wr.sf.95$gid),]
wr.sf.95$Runoff <- wr.sf.95$runoff *365*24*60*60 # from kg/m2/s to kg/m2/year or mm/y
wr.sf.95$logRunoff <- log(wr.sf.95$Runoff)
sfc_point_to_matrix = function(x) {
matrix(
unlist(x, use.names = FALSE),
nrow = length(x),
byrow = TRUE,
dimnames = list(1:length(x))
)
}
wr.centroid <- st_centroid(wr.sf.95, of.largest.polygon = T)
coord <- sfc_point_to_matrix(wr.centroid$geometry) %>% as.data.frame() %>% setNames(c("Longitude","Latitude"))
wr.sf.95 <- cbind(wr.sf.95, coord)
wr.sf.95$area <- st_area(wr.sf.95)
saveRDS(wr.sf.95, "WR/wr.sf.95.rds")noforest <- filter(wr.sf.95,forest == 0) #%>% dplyr::select(noforest,c("gid","area"))
#ggplot(noforest[which.max(noforest$area),])+geom_sf(aes(fill=as.factor(gid)))
weird.gid <- noforest$gid[which.max(noforest$area)]
weird.cat <- wr.sf.95 %>% filter(gid == weird.gid) %>% st_transform(CRS("EPSG:3035"))
corine <- raster::raster("Corine_Land_Cover/U2006_CLC2000_V2020_20u1.tif")
#weird.clc <- raster::extract(corine,weird.cat)
#saveRDS(weird.clc,"WR/clc/weird.clc.rds")
weird.clc <- readRDS("WR/clc/weird.clc.rds")
weird.tab <- sapply(weird.clc, function(x) tabulate(x,45))
weird.tabl.area <- weird.tab*prod(res(corine))
catchment.area <- colSums(weird.tabl.area)
weird.tab.prop <- sweep(weird.tabl.area,2,catchment.area, FUN = "/")
wr.sf.95$bogs[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[36]
wr.sf.95$forest[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[24]+weird.tab.prop[24]+weird.tab.prop[25]
wr.sf.95$arable[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[12]+weird.tab.prop[16]+weird.tab.prop[18]+weird.tab.prop[19]+weird.tab.prop[20]
wr.sf.95$glacier[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[34]
wr.sf.95$bare[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[31]+weird.tab.prop[32]+weird.tab.prop[33]
saveRDS(wr.sf.95,"WR/wr.sf.95.rds")wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
toreplace <- which(names(wr.sf.95) %in% c("ndvi","bogs","forest","arable","bare","water","temp","prec","alt","slope","tndep"))
names(wr.sf.95)[toreplace] <- c("NDVI","Bog","Forest","Arable","Bare","Water","Temp","Precip","Alt","Slope","TNdep")
wr.sf.95$area.x <- NULL
wr.sf.95$area.y <- NULL
saveRDS(wr.sf.95, "WR/wr.sf.95.rds")wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
m <- c("NDVI","Runoff","Precip","Temp","Bog","Forest","Arable","TNdep")
u <- c("NDVI on scale 0 to 1", "Mean Runoff, mm/y","Mean monthly precipitation,\n mm","Mean annual temperature,\n °C","Proportion of bogs","Proportion of Forest","Proportion of\narable land","Total Nitrogen Deposition (mg/m2)")
couleurs <- c("Greens","Blues","Blues","RdYlBu","YlOrBr","RdYlGn","RdYlGn","PuRd")
dir <- c(T,F,F,T,T,F,T,F)
for(i in m){
mysf <- wr.sf.95
map <- ggplot(mysf)+geom_sf(aes(fill = .data[[i]], col = .data[[i]]))+
scale_fill_continuous_divergingx(palette = "RdYlBu",
aesthetics = c("colour","fill"),
rev = dir[grep(i,m)],
name = paste(u[grep(i,m)],"\n(",i,")",sep=""))+
theme_minimal(base_size = 9)+
theme(legend.position = "bottom")
filename_i <- paste("WR/wr.map",i,"png",sep = ".")
ggsave(plot = map, filename = filename_i, dpi = "retina", width = 10, height = 12, units = "cm")
}
wr.sf.95 <- readRDS("WR/wr.sf.95.rds") %>% st_drop_geometry()
wr.sf.95$area_km2 <- wr.sf.95$area * 10^(-6)
wr.min <- lapply(wr.sf.95, min) %>% as.data.frame()
wr.max <- lapply(wr.sf.95, max) %>% as.data.frame()
wr.median <- lapply(wr.sf.95, median) %>% as.data.frame()
wr.mean <- lapply(wr.sf.95, mean) %>% as.data.frame()
wr.sd <- lapply(wr.sf.95, sd) %>% as.data.frame()
wr.sum <- rbind(wr.min, wr.max, wr.median, wr.mean, wr.sd)
rownames(wr.sum) <- c("Min","Max","Median","Mean","Sd")
wr.sum %>% dplyr::select(c("NDVI", "Forest", "Arable","Temp","Precip","TNdep","Runoff","area_km2")) %>% kable(digits = 2) %>%
kable_styling(bootstrap_options = "bordered")| NDVI | Forest | Arable | Temp | Precip | TNdep | Runoff | area_km2 | |
|---|---|---|---|---|---|---|---|---|
| Min | 0.11 | 0.00 | 0.00 | -4.24 | 434.55 | 63.88 | 55.17 | 27.3728 [m^2] |
| Max | 0.90 | 0.96 | 0.92 | 8.13 | 2842.75 | 2387.52 | 971.21 | 49832.2278 [m^2] |
| Median | 0.73 | 0.38 | 0.01 | 3.24 | 671.00 | 432.82 | 250.68 | 76.2258 [m^2] |
| Mean | 0.68 | 0.39 | 0.08 | 3.05 | 966.33 | 586.44 | 339.81 | 751.3046 [m^2] |
| Sd | 0.15 | 0.25 | 0.15 | 2.97 | 582.82 | 491.72 | 222.63 | 3752.6959 [m^2] |
2 Model prediction 1995
We model the mean TOC concentration in all the water region in 1995 based on the simple SEM model of the previous section, using NDVI, bogs and runoff as predictors.
2.1 Neighbour matrix
The neighbor matrix of the coastal drainage basins polygons were computed using their centroids.
wr.sf.95 <- readRDS("wr.sf.95.rds")
wr.centroid <- st_centroid(wr.sf.95, of_largest_polygon = T)
wr.neigh.set <- knearneigh(wr.centroid, k = 100)
wr.neigh.nb <- knn2nb(wr.neigh.set)
wr.kmat <- nb2listw(wr.neigh.nb)
saveRDS(wr.kmat,"WR/wr.kmat.rds")2.2 Model TOC
The TOC concentration in coastal drainage basins was modeled using the Spatial Error Linear Model with 5 predictors:
library(spatialreg)
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred <- predict(sem.fennoscandia.5,wr.sf.95,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
sem.pred.wr <- cbind(wr.sf.95,wr.pred$fit)
sem.pred.wr$TOC95 <- exp(sem.pred.wr$wr.pred.fit)
map.fit <- ggplot()+geom_sf(data = st_as_sf(sem.pred.wr), aes(fill = TOC95, col = TOC95))+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1,
name = "Modelled TOC 1995 (mg/L)")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
ggsave(plot = map.fit, filename = "WR/wr.toc.png",dpi = "retina", width = 6, height = 8, units = "cm")
map.log <- ggplot()+geom_sf(data = sem.pred.wr, aes(fill = wr.pred.fit, col = wr.pred.fit))+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1,
name = "Modelled log(TOC)")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
ggsave(plot = map.log, filename = "WR/wr.log.toc.png",dpi = "retina", width = 6, height = 8, units = "cm")
saveRDS(sem.pred.wr,"WR/sem.pred.wr.rds")knitr::include_graphics("WR/wr.log.toc.png")
Figure 6: Map of log(TOC) in coastal drainage basins for 1995
3 Forecast
3.1 Forecast of predictors
3.1.1 Forecast NDVI
3.1.1.1 Extract future temperature and precipitation
The forecast of the temperature and precipitation in the periods 2041-2060 and 2081-2100, for SSP 1-2.6 and SSP 3-7.0, were extracted based on the Future Worldclim dataset.
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.cnrm.temp <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp126_2041-2060.tif", band = 1)
ssp370.cnrm.temp <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp370_2041-2060.tif", band = 1)
ssp126.cnrm.temp.2050.crop <- raster::crop(ssp126.cnrm.temp,c(0,35,55,73))
png("worldclim_future/ssp126.temp.2050.png")
plot(ssp126.cnrm.temp.2050.crop)
dev.off()
ssp370.cnrm.temp.2050.crop <- raster::crop(ssp370.cnrm.temp,c(0,35,55,73))
png("worldclim_future/ssp370.temp.2050.png")
plot(ssp370.cnrm.temp.2050.crop)
dev.off()
ssp126.temp.wr <- raster::extract(ssp126.cnrm.temp,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm =T)
names(ssp126.temp.wr) <- c("gid", "temp126")
saveRDS(ssp126.temp.wr,"WR/ssp126.temp.wr.rds")
ssp370.temp.wr <- raster::extract(ssp370.cnrm.temp,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.temp.wr) <- c("gid", "temp370")
saveRDS(ssp370.temp.wr,"WR/ssp370.temp.wr.rds")wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.temp.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp126_2081-2100.tif", band = 1)
raster.ssp370.temp.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp370_2081-2100.tif", band = 1)
ssp126.cnrm.temp.2100.crop <- raster::crop(raster.ssp126.temp.2100,c(0,35,55,73))
png("worldclim_future/ssp126.temp.2100.png")
plot(ssp126.cnrm.temp.2100.crop)
dev.off()
ssp370.cnrm.temp.2100.crop <- raster::crop(raster.ssp370.temp.2100,c(0,35,55,73))
png("worldclim_future/ssp370.temp.2100.png")
plot(ssp370.cnrm.temp.2100.crop)
dev.off()
ssp126.temp.2100 <- raster::extract(raster.ssp126.temp.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm =T)
names(ssp126.temp.2100) <- c("gid", "temp126.2100")
saveRDS(ssp126.temp.2100,"WR/ssp126.temp.2100.rds")
ssp370.temp.2100 <- raster::extract(raster.ssp370.temp.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.temp.2100) <- c("gid", "temp370.2100")
saveRDS(ssp370.temp.2100,"WR/ssp370.temp.2100.rds")listgrob <- list("worldclim_future/ssp126.temp.2050.png", "worldclim_future/ssp370.temp.2050.png","worldclim_future/ssp126.temp.2100.png","worldclim_future/ssp370.temp.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("a) SSP126-2050","b) SSP370-2050","c) SSP126-2100","d) SSP370-2100"), label_size = 10)
Figure 7: Map of forecasted temperature (ºC) for 2050 and 2100, under SSP 1-2.6 and SSP 3-7.0
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.cnrm.prec <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp126_2041-2060.tif", band = 12)
ssp370.cnrm.prec <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp370_2041-2060.tif", band = 12)
ssp126.cnrm.prec.crop <- raster::crop(ssp126.cnrm.prec,c(0,35,55,73))
png("worldclim_future/ssp126.2050.prec.png")
plot(ssp126.cnrm.prec.crop)
dev.off()
ssp370.cnrm.prec.crop <- raster::crop(ssp370.cnrm.prec,c(0,35,55,73))
png("worldclim_future/ssp370.2050.prec.png")
plot(ssp370.cnrm.prec.crop)
dev.off()
ssp126.prec.wr <- raster::extract(ssp126.cnrm.prec,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp126.prec.wr) <- c("gid", "prec126")
saveRDS(ssp126.prec.wr,"WR/ssp126.prec.wr.rds")
ssp370.prec.wr <- raster::extract(ssp370.cnrm.prec,dplyr::select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.prec.wr) <- c("gid", "prec370")
saveRDS(ssp370.prec.wr,"WR/ssp370.prec.wr.rds")wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.prec.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp126_2081-2100.tif", band = 12)
raster.ssp370.prec.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp370_2081-2100.tif", band = 12)
ssp126.cnrm.prec.2100.crop <- raster::crop(raster.ssp126.prec.2100,c(0,35,55,73))
png("worldclim_future/ssp126.2100.prec.png")
plot(ssp126.cnrm.prec.2100.crop)
dev.off()
png("worldclim_future/ssp370.2100.prec.png")
ssp370.cnrm.prec.2100.crop <- raster::crop(raster.ssp370.prec.2100,c(0,35,55,73))
plot(ssp370.cnrm.prec.2100.crop)
ssp126.prec.2100 <- raster::extract(raster.ssp126.prec.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp126.prec.2100) <- c("gid", "prec126.2100")
saveRDS(ssp126.prec.2100,"WR/ssp126.prec.2100.rds")
ssp370.prec.2100 <- raster::extract(raster.ssp370.prec.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.prec.2100) <- c("gid", "prec370.2100")
saveRDS(ssp370.prec.2100,"WR/ssp370.prec.2100.rds")wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.temp.wr <- readRDS("WR/ssp126.temp.wr.rds")
ssp370.temp.wr <- readRDS("WR/ssp370.temp.wr.rds")
ssp126.prec.wr <- readRDS("WR/ssp126.prec.wr.rds")
ssp370.prec.wr <- readRDS("WR/ssp370.prec.wr.rds")
ssp126.temp.2100 <- readRDS("WR/ssp126.temp.2100.rds")
ssp370.temp.2100 <- readRDS("WR/ssp370.temp.2100.rds")
ssp126.prec.2100 <- readRDS("WR/ssp126.prec.2100.rds")
ssp370.prec.2100 <- readRDS("WR/ssp370.prec.2100.rds")
# SSP126
df.ssp126 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp126.temp.wr, by = "gid") %>% merge(ssp126.prec.wr, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp126$logPrecip <- df.ssp126$Precip %>% log()
saveRDS(df.ssp126, "WR/df.ssp126.rds")
# SSP 370
df.ssp370 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp370.temp.wr, by = "gid") %>% merge(ssp370.prec.wr, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp370$logPrecip <- df.ssp370$Precip %>% log()
saveRDS(df.ssp370, "WR/df.ssp370.rds")
# SSP126 2100
df.ssp126.2100 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp126.temp.2100, by = "gid") %>% merge(ssp126.prec.2100, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp126.2100$logPrecip <- df.ssp126.2100$Precip %>% log()
saveRDS(df.ssp126.2100, "WR/df.ssp126.2100.rds")
# SSP 370 2100
df.ssp370.2100 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp370.temp.2100, by = "gid") %>% merge(ssp370.prec.2100, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp370.2100$logPrecip <- df.ssp370.2100$Precip %>% log()
saveRDS(df.ssp370.2100, "WR/df.ssp370.2100.rds")
listgrob <- list("worldclim_future/ssp126.2050.prec.png", "worldclim_future/ssp370.2050.prec.png","worldclim_future/ssp126.2100.prec.png","worldclim_future/ssp370.2100.prec.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("a) SSP126-2050","b) SSP370-2050","c) SSP126-2100","d) SSP370-2100"), label_size = 10)
Figure 8: Map of forecasted precipitations (mm/s) for 2050 and 2100, under SSP 1-2.6 and SSP 3-7.0
3.1.1.2 NDVI model based on past temperature and precipitation values
We model NDVI based on temperature and precipitation, and use this model to predict the future NDVI in 2050. Choice of future data: ssp 126 as best scenario, ssp 370 as worse scenario (Hausfather 2019; CORDEX 2021). CRNM according to Raju and Kumar (2020).
Using a SELM model for NDVI, we notice that high residuals are observed in catchments where original NDVI is lower than 0,4. Myers-Smith et al. (2020) show that the relationship between NDVI and the Arctic biomass is not linear: under 0,4, an increase in NDVI reflects the vegetalization of bare ground, while over 0,4, an increase of NDVI reflects an increased productivity in established forests.
Several models were tested to model the NDVI: a linear model, a spatial error linear model, a polynomial model, and a polynomial model combined with a binomial model. The NDVI was fitted on the 1995 data from the Northern Lakes Survey and tested on on the coastal drainage basins.
library(betareg)
library(mgcv)
library(VGAM)
wr.kmat <- readRDS("WR/wr.kmat.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
fennoscandia <- readRDS("fennoscandia.33.rds")
fennoscandia$Latitude <- fennoscandia$latitude
fen.kmat <- readRDS("k.neigh.w.rds")
ndvi.cor <- cor(dplyr::select(st_drop_geometry(fennoscandia), c("NDVI", "Temp", "Precip", "Arable", "Bare", "Latitude", "Slope","Runoff")))
corrplot::corrplot(ndvi.cor, method = "number")
par(mfrow = c(2,3))
plot(NDVI~Temp+log(Precip)+Arable+Bare+Latitude+Slope, data = st_drop_geometry(fennoscandia))
par(mfrow = c(3,2))
#LM
summary(m1 <- lm(NDVI ~ Temp + log(Precip), data = fennoscandia))
summary(m1 <- step(m1))
plot(predict(m1), fennoscandia$NDVI, pch = ".", main="a ) NDVI - LM", ylab="Observed NDVI", xlab="Predicted",abline(0,1))
text(x = 0.8, y = 0.4, label = paste("r2 =", round(summary(m1)$adj.r.squared,2), sep= ""))
#SELM
fen.kmat <- readRDS("k.neigh.w.rds")
selm.model <- errorsarlm(formula = NDVI~Temp+logPrecip, data = fennoscandia, listw = fen.kmat)
saveRDS(selm.model, "WR/selm.ndvi.rds")
selm.model <- readRDS("WR/selm.ndvi.rds")
AIC <- 2*(3-1)-2*ndvi.model$LL
plot(ndvi.model$fitted.value, fennoscandia$NDVI, pch = ".", main="b) NDVI - SELM", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.8, y = 0.5, label = paste("r2 =", round(cor(ndvi.model$fitted.values, fennoscandia$NDVI)^2,2), sep= ""))
# quadratic
summary(m2 <- lm(NDVI ~ (Temp + I(Temp^2)) * (log(Precip) + I(log(Precip)^2)), data = fennoscandia))
summary(m2 <- step(m2))
plot(predict(m2), fennoscandia$NDVI, pch = ".", main="c) NDVI - quadratic model", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.7, y = 0.4, label = paste("r2 =", round(summary(m2)$adj.r.squared,2), sep= ""))
# poly
summary(m3 <- lm(NDVI ~ poly(cbind(Temp, logPrecip), degree = 2), data = fennoscandia))
summary(m3 <- step(m3))
plot(predict(m3), fennoscandia$NDVI, pch = ".", main="d) NDVI - polynomial ", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.7, y = 0.4, label = paste("r2 =", round(summary(m3)$adj.r.squared,2), sep= ""))
saveRDS(m3,"WR/poly.ndvi.rds")
# glm
summary(m4 <- glm(NDVI ~ poly(cbind(Temp, logPrecip), degree = 2), data = fennoscandia, family = "quasibinomial"))
#summary(m4 <- step(m4))
plot(predict(m4, type = "response"), fennoscandia$NDVI, pch = ".", main="e) NDVI - glm ", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.7, y = 0.4, label = paste("r =", cor(fennoscandia$NDVI, m4$fitted.values), sep= ""))
saveRDS(m4, "WR/glm.ndvi.rds")
#betareg
summary(m5 <- betareg(NDVI ~ poly(cbind(Temp, logPrecip), degree = 2), data = fennoscandia, link = "logit"))
plot(predict(m5, type = "response"), fennoscandia$NDVI, pch = ".", main="e) NDVI - glm ", ylab="Observed", xlab="Predicted",abline(0,1))
saveRDS(m5, "WR/betareg.ndvi.rds")
Figure 9: Test of 4 NDVI models, a) linear model, b) spatial error linear model, c) quadratic model, d) polynomial model
The SELM and the polynomial model had the best results. We tested them on the coastal drainage basins dataset.
k.neigh.w <- readRDS("k.neigh.w.rds")
fennoscandia <- readRDS("fennoscandia.rds")
fennoscandia.33 <- readRDS("fennoscandia.33.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
selm.ndvi <- readRDS("WR/selm.ndvi.rds")
selm.ndvi.fen <- predict(selm.ndvi,wr.sf.95,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.selm.ndvi.model <- wr.sf.95 %>% st_as_sf() %>% mutate(ndvi.fit = selm.ndvi.fen$fit, ndvi.residuals = wr.sf.95$NDVI - selm.ndvi.fen$fit)
wr.selm.ndvi.model$area <- st_area(wr.selm.ndvi.model) %>% as.numeric()
ggplot(wr.selm.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.selm.ndvi.model$NDVI,wr.selm.ndvi.model$ndvi.fit)^2,2)), x = 0.25, y = 0.5, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for SELM")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.selm.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals.selm <- ggplot()+
geom_sf(data = wr.selm.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "SELM residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for SELM")
ggsave(plot = wr.residuals.selm, filename = "WR/wr.selm.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")poly.ndvi <- readRDS("WR/poly.ndvi.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
fennoscandia <- readRDS("fennoscandia.33.rds")
poly.ndvi.wr <- predict(poly.ndvi,wr.sf.95, type = "response")
wr.poly.ndvi.model <- wr.sf.95 %>% mutate(ndvi.fit = poly.ndvi.wr, ndvi.residuals = wr.sf.95$NDVI - poly.ndvi.wr)
ggplot(wr.poly.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.poly.ndvi.model$NDVI,wr.poly.ndvi.model$ndvi.fit),2)), x = 0.75, y = 0.4, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for polynomial model")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.poly.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals <- ggplot()+
geom_sf(data = wr.poly.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "Polynomial residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for Polynomial model")
ggsave(plot = wr.residuals, filename = "WR/wr.poly.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")
glm.ndvi <- readRDS("WR/glm.ndvi.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
glm.ndvi.wr <- predict(glm.ndvi,wr.sf.95, type = "response")
wr.glm.ndvi.model <- wr.sf.95 %>% mutate(ndvi.fit = glm.ndvi.wr, ndvi.residuals = wr.sf.95$NDVI - glm.ndvi.wr)
ggplot(wr.glm.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.glm.ndvi.model$NDVI,wr.glm.ndvi.model$ndvi.fit),2)), x = 0.8, y = 0.3, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for binomial model")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.glm.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals <- ggplot()+
geom_sf(data = wr.glm.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "glm residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for Binomial model")
ggsave(plot = wr.residuals, filename = "WR/wr.glm.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")
betareg.ndvi <- readRDS("WR/betareg.ndvi.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
betareg.ndvi.wr <- predict(betareg.ndvi,wr.sf.95, type = "response")
wr.betareg.ndvi.model <- wr.sf.95 %>% mutate(ndvi.fit = betareg.ndvi.wr, ndvi.residuals = wr.sf.95$NDVI - betareg.ndvi.wr)
ggplot(wr.betareg.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.betareg.ndvi.model$NDVI,wr.betareg.ndvi.model$ndvi.fit),2)), x = 0.8, y = 0.3, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for binomial model")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.betareg.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals <- ggplot()+
geom_sf(data = wr.betareg.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "betareg residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for Binomial model")
ggsave(plot = wr.residuals, filename = "WR/wr.betareg.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")listgrob1 <- list("WR/wr.ndvi.selm.results.png","WR/wr.ndvi.poly.results.png","WR/wr.ndvi.glm.results.png","WR/wr.ndvi.betareg.results.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob1, ncol = 2, nrow = 2, labels = c("a)","b)","c)","d)"))
listgrob2 <- list("WR/wr.selm.ndvi.residuals.png","WR/wr.poly.ndvi.residuals.png", "WR/wr.glm.ndvi.residuals.png","WR/wr.betareg.ndvi.residuals.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob2, ncol= 2, nrow = 2,labels = c("e)","f)","g)","h)"))
Figure 10: Results of NDVI models (SELM and polynomial): observed vs predicted for a) SELM and b) polynoimal model, and map of residuals (observed - predicted) for c) SELM and d) Polynomial
selm.ndvi <- readRDS("WR/selm.ndvi.rds")
poly.ndvi <- readRDS("WR/poly.ndvi.rds")
glm.ndvi <- readRDS("WR/glm.ndvi.rds")
betareg.ndvi <- readRDS("WR/betareg.ndvi.rds")
fennoscandia <- readRDS("fennoscandia.33.rds")
xtemp <- seq(min(fennoscandia$Temp), max(fennoscandia$Temp), 0.1)
yprecip <- seq(min(fennoscandia$logPrecip), max(fennoscandia$logPrecip), 0.1)
data.fit <- expand.grid(Temp = xtemp, logPrecip = yprecip)
mtrx3d <- predict(selm.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp1 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot SELM")+
theme_minimal(base_size = 10)
mtrx3d <- predict(poly.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp2 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot polynomial model")+
theme_minimal(base_size = 10)
mtrx3d <- predict(glm.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp3 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot binomial model")+
theme_minimal(base_size = 10)
mtrx3d <- predict(betareg.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp4 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..), col = "black")+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot beta model")+
theme_minimal(base_size = 10)
grid.arrange(cp1,cp2,cp3,cp4, nrow = 2, ncol = 2)
The binomial model preforms better: low residuals and prediction range between 0 and 1. This model will be used to forecast NDVI in 2050 and in 2100.
3.1.1.3 Forecast NDVI in 2050 and 2100
df.ssp126 <- readRDS("WR/df.ssp126.rds")
df.ssp370 <- readRDS("WR/df.ssp370.rds")
df.ssp126.2100 <- readRDS("WR/df.ssp126.2100.rds")
df.ssp370.2100 <- readRDS("WR/df.ssp370.2100.rds")
ndvi.model <- readRDS("WR/betareg.ndvi.rds")
wr.kmat <- readRDS("WR/wr.kmat.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- wr.sf.95$Precip %>% log()
#NDVI 126 2050
# ndvi.ssp126.fit <- predict(object = ndvi.model, newdata = df.ssp126, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
# ndvi.ssp126 <- cbind(df.ssp126, ndvi.ssp126.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp126.fit <- predict(object = ndvi.model, newdata = df.ssp126, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp126 <- cbind(df.ssp126, ndvi.ssp126.fit)
saveRDS(ndvi.ssp126,"WR/ndvi.ssp126.rds")
# NDVI 370 2050
# ndvi.ssp370.fit <- predict(object = ndvi.model, newdata = df.ssp370, listw = wr.kmat, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
#
# ndvi.ssp370 <- cbind(df.ssp370, ndvi.ssp370.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp370.fit <- predict(object = ndvi.model, newdata = df.ssp370, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp370 <- cbind(df.ssp370, ndvi.ssp370.fit)
saveRDS(ndvi.ssp370,"WR/ndvi.ssp370.rds")
#NDVI 126 2010
# ndvi.ssp126.2100.fit <- predict(object = ndvi.model, newdata = df.ssp126.2100, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
#
# ndvi.ssp126.2100 <- cbind(df.ssp126.2100, ndvi.ssp126.2100.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp126.2100.fit <- predict(object = ndvi.model, newdata = df.ssp126.2100, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp126.2100 <- cbind(df.ssp126.2100, ndvi.ssp126.2100.fit)
saveRDS(ndvi.ssp126.2100,"WR/ndvi.ssp126.2100.rds")
#NDVI 370 2100
# ndvi.ssp370.2100.fit <- predict(object = ndvi.model, newdata = df.ssp370.2100, listw = wr.kmat, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% setNames(c("ndvi","trend","signal"))
# ndvi.ssp370.2100 <- cbind(df.ssp370.2100, ndvi.ssp370.2100.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp370.2100.fit <- predict(object = ndvi.model, newdata = df.ssp370.2100, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp370.2100 <- cbind(df.ssp370.2100, ndvi.ssp370.2100.fit)
saveRDS(ndvi.ssp370.2100,"WR/ndvi.ssp370.2100.rds")
Figure 10: Results of NDVI models (SELM and polynomial): observed vs predicted for a) SELM and b) polynoimal model, and map of residuals (observed - predicted) for c) SELM and d) Polynomial
3.1.2 Forecast runoff 2050 and 2100
In the historical dataset in Worldclim, the runoff is given as a yearly average in kg/m2/s, or mm/s.
In the CNRM-CM6 forecast we get the monthly forecast for runoff (in mm/s) and convert it in mm/y.
listgrob <- list("CORDEX_future/runoff.ssp126.2050.png", "CORDEX_future/runoff.ssp370.2050.png","CORDEX_future/runoff.ssp126.2100.png","CORDEX_future/runoff.ssp370.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
ggarrange(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("SSP126-2050","SSP370-2050","SSP126-2100","SSP370-2100"))
wr.sf <- readRDS("WR/wr.sf.rds")
ssp126.runoff.wr <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp126_r1i1p1f2_gr_201501-210012.nc", bands = c(313:552))
ssp126.runoff.wr.mean <- ssp126.runoff.wr %>% mean()
ssp126.runoff.wr.mean.crop <- raster::crop(ssp126.runoff.wr.mean,c(0,35,55,73))
plot(ssp126.runoff.wr.mean.crop)
ssp126.runoff <- raster::extract(ssp126.runoff.wr.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp126.runoff) <- c("gid","runoff")
ssp126.runoff$Runoff <- ssp126.runoff$runoff *365*24*60*60
saveRDS(ssp126.runoff,"WR/ssp126.runoff.2050.rds")
ssp370.runoff.wr <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp370_r1i1p1f2_gr_201501-210012.nc", bands = c(313:552))
ssp370.runoff.wr.mean <- ssp370.runoff.wr %>% mean()
ssp370.runoff.wr.mean.crop <- raster::crop(ssp370.runoff.wr.mean,c(0,35,55,73))
plot(ssp370.runoff.wr.mean.crop)
ssp370.runoff <- raster::extract(ssp370.runoff.wr.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp370.runoff) <- c("gid","runoff")
ssp370.runoff$Runoff <- ssp370.runoff$runoff *365*24*60*60
saveRDS(ssp370.runoff,"WR/ssp370.runoff.2050.rds")wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.runoff.2100 <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp126_r1i1p1f2_gr_201501-210012.nc", bands = c(805:1032))
raster.ssp126.runoff.2100.mean <- raster.ssp126.runoff.2100 %>% mean()
raster.ssp126.runoff.2100.mean.crop <- raster::crop(raster.ssp126.runoff.2100.mean,c(0,35,55,73))
plot(raster.ssp126.runoff.2100.mean.crop)
ssp126.runoff.2100 <- raster::extract(raster.ssp126.runoff.2100.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp126.runoff.2100) <- c("gid","runoff")
ssp126.runoff.2100$Runoff <- ssp126.runoff.2100$runoff *365*24*60*60
saveRDS(ssp126.runoff.2100,"WR/ssp126.runoff.2100.rds")
raster.ssp370.runoff.2100 <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp370_r1i1p1f2_gr_201501-210012.nc", bands = c(805:1032))
raster.ssp370.runoff.2100.mean <- raster.ssp370.runoff.2100 %>% mean()
raster.ssp370.runoff.2100.mean.crop <- raster::crop(raster.ssp370.runoff.2100.mean,c(0,35,55,73))
plot(raster.ssp370.runoff.2100.mean.crop)
ssp370.runoff.2100 <- raster::extract(raster.ssp370.runoff.2100.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp370.runoff.2100) <- c("gid","runoff")
ssp370.runoff.2100$Runoff <- ssp370.runoff.2100$runoff *365*24*60*60
saveRDS(ssp370.runoff.2100,"WR/ssp370.runoff.2100.rds")3.1.3 Future TNdep
Future TNdep will largely depend on the public policies implemented under the scenarios SSP 1-2.6 and 3-7.0.
TNdep data are from 2000.
Emission reduction between 2005 and 2020 (European Environment Agency 2021):
- NH3: -8%
- NOx -42%
Predictions for future emissions for NOx (summary for policymakers AR6 (Masson-Delmotte et al. 2021)):
- SSP1-2.6: from 11 MtN2O/y to 9 MtN2O (18% decrease) in 2050, stable in 2100 (also 20% compared to 2015)
- SSP3-7.0: from 11 MtN2o/y to 15 MtN2O (137% increase) in 2050, and 21 MtN2O (190% increase) in 2100
Assumption for NH3:
- SSP1-2.6: Decrease according to targets (-20%)
- SSP3-7.0: No decrease
ndep.wr <- readRDS("WR/ndep.wr.rds")
ndep.sum <- ndep.wr@data %>% apply(2, sum)
(ndep.sum["woxn"]+ndep.sum["doxn"])/ndep.sum["tndep"]
(ndep.sum["wrdn"]+ndep.sum["drdn"])/ndep.sum["tndep"]
# SSP126 2050
ssp126.ndep.2050 <- data.frame(gid = ndep.wr$gid)
ssp126.ndep.2050$oxn <- ndep.wr$woxn*0.55*0.8 + ndep.wr$doxn*0.55*0.8 # -45% for decrease between 2000 and 2020, -20% between 2020 and 2050
ssp126.ndep.2050$rdn <- ndep.wr$wrdn*0.9*0.8 + ndep.wr$drdn*0.9*0.8 # -10% for decrease between 2000 and 2020, -20% between 2020 and 2050
ssp126.ndep.2050$TNdep <- ssp126.ndep.2050$oxn + ssp126.ndep.2050$rdn
saveRDS(ssp126.ndep.2050,"WR/ssp126.ndep.2050.rds")
# SSP126 2100
ssp126.ndep.2100 <- data.frame(gid = ndep.wr$gid)
ssp126.ndep.2100$oxn <- ndep.wr$woxn*0.55*0.8+ndep.wr$doxn*0.55*0.8 # -45% for decrease between 2000 and 2020, -20% between 2020 and 2100
ssp126.ndep.2100$rdn <- ndep.wr$wrdn*0.9*0.8+ndep.wr$drdn*0.9*0.8 # -10% for decrease between 2000 and 2020, -20% between 2020 and 2100
ssp126.ndep.2100$TNdep <- ssp126.ndep.2100$oxn + ssp126.ndep.2100$rdn
saveRDS(ssp126.ndep.2100,"WR/ssp126.ndep.2100.rds")
# SSP370 2050
ssp370.ndep.2050 <- data.frame(gid = ndep.wr$gid)
ssp370.ndep.2050$oxn <- ndep.wr$woxn*0.55*1.37+ndep.wr$doxn*0.55*1.37 # -45% for decrease between 2000 and 2020, 137% increase between 2020 and 2100
ssp370.ndep.2050$rdn <- ndep.wr$wrdn*0.9+ndep.wr$drdn*0.9 # -10% for decrease between 2000 and 2020, no change between 2020 and 2100
ssp370.ndep.2050$TNdep <- ssp370.ndep.2050$oxn + ssp370.ndep.2050$rdn
saveRDS(ssp370.ndep.2050,"WR/ssp370.ndep.2050.rds")
# SSP370 2100
ssp370.ndep.2100 <- data.frame(gid = ndep.wr$gid)
ssp370.ndep.2100$oxn <- ndep.wr$woxn*0.55*1.9+ndep.wr$doxn*0.55*1.9 # -45% for decrease between 2000 and 2020, 190% increase between 2020 and 2100
ssp370.ndep.2100$rdn <- ndep.wr$wrdn*0.9+ndep.wr$drdn*0.9 # -10% for decrease between 2000 and 2020, no change between 2020 and 2100
ssp370.ndep.2100$TNdep <- ssp370.ndep.2100$oxn + ssp370.ndep.2100$rdn
saveRDS(ssp370.ndep.2100,"WR/ssp370.ndep.2100.rds")3.1.4 Gather datasets and plot
The datasets were gathered into a single dataframe.
ndvi.ssp126 <- readRDS("WR/ndvi.ssp126.rds")
ndvi.ssp370 <- readRDS("WR/ndvi.ssp370.rds")
ssp126.runoff.2050 <- readRDS("WR/ssp126.runoff.2050.rds")
ssp370.runoff.2050 <- readRDS("WR/ssp370.runoff.2050.rds")
ssp126.ndep.2050 <- readRDS("WR/ssp126.ndep.2050.rds")
ssp370.ndep.2050 <- readRDS("WR/ssp370.ndep.2050.rds")
ndvi.ssp126.2100 <- readRDS("WR/ndvi.ssp126.2100.rds")
ndvi.ssp370.2100 <- readRDS("WR/ndvi.ssp370.2100.rds")
ssp126.runoff.2100 <- readRDS("WR/ssp126.runoff.2100.rds")
ssp370.runoff.2100 <- readRDS("WR/ssp370.runoff.2100.rds")
ssp126.ndep.2100 <- readRDS("WR/ssp126.ndep.2100.rds")
ssp370.ndep.2100 <- readRDS("WR/ssp370.ndep.2100.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
# 2050 -----
wr.sf.126.2050 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp126, !c("Arable", "Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp126.runoff.2050)), by = "gid") %>% merge(ssp126.ndep.2050)
wr.sf.370.2050 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp370, !c("Arable","Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp370.runoff.2050)), by = "gid") %>% merge(ssp370.ndep.2050)
na.lines.126 <- which(complete.cases(st_drop_geometry(wr.sf.126.2050))==F)
na.lines.370 <- which(complete.cases(st_drop_geometry(wr.sf.370.2050))==F)
wr.sf.126.2050 <- wr.sf.126.2050[-na.lines.126,]
wr.sf.370.2050 <- wr.sf.370.2050[-na.lines.370,]
wr.sf.126.2050$logRunoff <- wr.sf.126.2050$Runoff %>% log()
wr.sf.370.2050$logRunoff <- wr.sf.370.2050$Runoff %>% log()
saveRDS(wr.sf.126.2050, "WR/wr.sf.126.2050.rds")
saveRDS(wr.sf.370.2050, "WR/wr.sf.370.2050.rds")
# 2100 ----
wr.sf.126.2100 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp126.2100, !c("Arable","Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp126.runoff.2100)), by = "gid") %>% merge(ssp126.ndep.2100)
wr.sf.370.2100 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp370.2100, !c("Arable","Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp370.runoff.2100)), by = "gid") %>% merge(ssp370.ndep.2100)
na.lines.126 <- which(complete.cases(st_drop_geometry(wr.sf.126.2100))==F)
na.lines.370 <- which(complete.cases(st_drop_geometry(wr.sf.370.2100))==F)
wr.sf.126.2100 <- wr.sf.126.2100[-na.lines.126,]
wr.sf.370.2100 <- wr.sf.370.2100[-na.lines.370,]
wr.sf.126.2100$logRunoff <- wr.sf.126.2100$Runoff %>% log()
wr.sf.370.2100$logRunoff <- wr.sf.370.2100$Runoff %>% log()
saveRDS(wr.sf.126.2100, "WR/wr.sf.126.2100.rds")
saveRDS(wr.sf.370.2100, "WR/wr.sf.370.2100.rds")
Then we plot the difference between the forecasted variables in 2050 and 2100, compared to 1995.
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
t1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 1-2.6, 2050", x = "Temperature 1995 (°C)")+theme_minimal()
t2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 1-2.6, 2100", x = "Temperature 1995 (°C)")+theme_minimal()
t3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 3-7.0, 2050", x = "Temperature 1995 (°C)")+theme_minimal()
t4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 3-7.0, 2100", x = "Temperature 1995 (°C)")+theme_minimal()
ggarrange(t1,t2,t3,t4, nrow = 2, ncol = 2)
p1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Precipitations SSP 1-2.6, 2050", x = "Precipitations 1995 (mm/y)")+theme_minimal()
p2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y= "Precipitations SSP 1-2.6, 2100", x = "Precipitations 1995 (mm/y)")+theme_minimal()
p3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Precipitations SSP 3-7.0, 2050", x = "Precipitations 1995 (mm/y)")+theme_minimal()
p4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Precipitations SSP 3-7.0, 2100", x = "Precipitations 1995 (mm/y)")+theme_minimal()
ggarrange(p1,p2,p3,p4, nrow = 2, ncol = 2)
# d1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
# d2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
# d3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
# d4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
#
# ggarrange(d1,d2,d3,d4, nrow = 2, ncol = 2)
Figure 11: Forecasted temperature and precipitation compared to 1995 values
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
n1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = NDVI95, y = NDVI))+ geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 1-2.6, 2050")+ theme_minimal()
n2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = NDVI95, y = NDVI))+geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 1-2.6, 2100")+theme_minimal()
n3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = NDVI95, y = NDVI))+geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 3-7.0, 2050")+theme_minimal()
n4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = NDVI95, y = NDVI))+geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 3-7.0, 2100")+theme_minimal()
ggarrange(n1,n2,n3,n4, nrow = 2, ncol = 2)
Figure 11: Forecasted temperature and precipitation compared to 1995 values
glm.ndvi <- readRDS("WR/glm.ndvi.rds")
fennoscandia <- readRDS("fennoscandia.33.rds")
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
par(mfrow = c(2,2))
xtemp1 <- seq(min(wr.sf.126.2050$Temp), max(wr.sf.126.2050$Temp), 0.1)
yprecip1 <- seq(min(wr.sf.126.2050$logPrecip), max(wr.sf.126.2050$logPrecip), 0.1)
data.fit1 <- expand.grid(Temp = xtemp1, logPrecip = yprecip1)
mtrx3d1 <- predict(glm.ndvi, newdata = data.fit1, type = "response")
mtrx.melt1 <- cbind(data.fit1, mtrx3d1)
names(mtrx.melt1) <- c("Temp","logPrecip","NDVI")
fcp1 <- ggplot(mtrx.melt1,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.126.2050, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 1-2.6, 2050")
xtemp2 <- seq(min(wr.sf.126.2100$Temp), max(wr.sf.126.2100$Temp), 0.1)
yprecip2 <- seq(min(wr.sf.126.2100$logPrecip), max(wr.sf.126.2100$logPrecip), 0.1)
data.fit2 <- expand.grid(Temp = xtemp2, logPrecip = yprecip2)
mtrx3d2 <- predict(glm.ndvi, newdata = data.fit2, type = "response")
mtrx.melt2 <- cbind(data.fit2, mtrx3d2)
names(mtrx.melt2) <- c("Temp","logPrecip","NDVI")
fcp2 <- ggplot(mtrx.melt2,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.126.2100, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 1-2.6, 2100")
xtemp3 <- seq(min(wr.sf.370.2050$Temp), max(wr.sf.370.2050$Temp), 0.1)
yprecip3 <- seq(min(wr.sf.370.2050$logPrecip), max(wr.sf.370.2050$logPrecip), 0.1)
data.fit3 <- expand.grid(Temp = xtemp3, logPrecip = yprecip3)
mtrx3d3 <- predict(glm.ndvi, newdata = data.fit3, type = "response")
mtrx.melt3 <- cbind(data.fit3, mtrx3d3)
names(mtrx.melt3) <- c("Temp","logPrecip","NDVI")
fcp3 <- ggplot(mtrx.melt3,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.370.2050, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 3-7.0, 2050")
xtemp4 <- seq(min(wr.sf.370.2100$Temp), max(wr.sf.370.2100$Temp), 0.1)
yprecip4 <- seq(min(wr.sf.370.2100$logPrecip), max(wr.sf.370.2100$logPrecip), 0.1)
data.fit4 <- expand.grid(Temp = xtemp4, logPrecip = yprecip4)
mtrx3d4 <- predict(glm.ndvi, newdata = data.fit4, type = "response")
mtrx.melt4 <- cbind(data.fit4, mtrx3d4)
names(mtrx.melt4) <- c("Temp","logPrecip","NDVI")
fcp4 <- ggplot(mtrx.melt4,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.370.2100, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 3-7.0, 2100")
grid.arrange(fcp1,fcp2,fcp3,fcp4,nrow=2,ncol=2)
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
mn0 <- ggplot(wr.sf.126.2050)+geom_sf(aes(fill = NDVI95, col = NDVI95)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI 1995", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn0, filename = "WR/NDVI_1995.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn1 <- ggplot(wr.sf.126.2050)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP126-2050", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn1, filename = "WR/NDVI_SSP126_2050.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn2 <- ggplot(wr.sf.126.2100)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP126-2100", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn2, filename = "WR/NDVI_SSP126_2100.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn3 <- ggplot(wr.sf.370.2050)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP370-2050", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn3, filename = "WR/NDVI_SSP370_2050.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn4 <- ggplot(wr.sf.370.2100)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP370-2100", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn4, filename = "WR/NDVI_SSP370_2100.png",dpi = "retina", width = 12, height = 16, units = "cm")
listgrob <- list("WR/NDVI_1995.png","white.png","WR/NDVI_SSP126_2050.png", "WR/NDVI_SSP126_2100.png","WR/NDVI_SSP370_2050.png","WR/NDVI_SSP370_2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 3, align = "v", greedy = T)
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
diff.var <- function(sf, vars1, vars2){
sf2 <- sf
for (i in 1:length(vars2)){
newname <- paste("diff", vars2[i], sep = "")
x <- sf2[,vars2[i]] %>% st_drop_geometry()
y <- sf2[,vars1[i]] %>% st_drop_geometry()
z <- x-y
names(z) <- newname
sf2 <- cbind(sf2, z)
}
return(sf2)
}
wr.sf.126.2050 <- diff.var(wr.sf.126.2050, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.370.2050 <- diff.var(wr.sf.370.2050, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.126.2100 <- diff.var(wr.sf.126.2100, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.370.2100 <- diff.var(wr.sf.370.2100, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
list.sf <- list(wr.sf.126.2050, wr.sf.126.2100, wr.sf.370.2050, wr.sf.370.2100)
names(list.sf) <- c("SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
all.diff <- data.frame(gid = c(wr.sf.126.2050$gid, wr.sf.370.2050$gid, wr.sf.126.2100$gid, wr.sf.370.2100$gid))
for(x in c("diffTemp","diffPrecip","diffNDVI","diffRunoff","diffTNdep")){
all.vec <- c()
for(i in 1:length(list.sf)){
vec <- select(st_drop_geometry(list.sf[[i]]),all_of(x))
all.vec <<- rbind(all.vec,vec)
}
all.diff <- cbind(all.diff, all.vec)
}
max.lim <- apply(all.diff, 2, max, na.rm = T)
min.lim <- apply(all.diff, 2, min, na.rm = T)
# x1 <- filter(wr.sf.126.2050, NDVI95 > 0.4)
# x2 <- filter(wr.sf.370.2050, NDVI95 > 0.4)
# x3 <- filter(wr.sf.126.2100, NDVI95 > 0.4)
# x4 <- filter(wr.sf.370.2100, NDVI95 > 0.4)
# all.ndvi <- c(x1$NDVI/x1$NDVI95*100-100,x2$NDVI/x2$NDVI95*100-100, x3$NDVI/x3$NDVI95*100-100, x4$NDVI/x4$NDVI95*100-100)
# palettes <- data.frame(diffTemp = brewer.pal(name = "RdYlBu", n = 9)[c(1,5,9)],
# diffPrecip = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)],
# diffNDVI = brewer.pal(name = "RdYlGn", n = 9)[c(9,5,1)],
# diffRunoff = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)],
# diffTNdep = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)])
legends <- c("Difference Temperature (°C)", "Difference Precipitation (mm/y)", "Difference NDVI", "Difference Runoff (mm/y)","DIfference Nitrogen Deposition (ug/m2)")
diffs <- c("diffTemp","diffPrecip","diffNDVI","diffRunoff","diffTNdep")
rev.pal <- c(T,F,T,F,T)
for(i in c(1:length(list.sf))){
for(x in diffs[3]){
g <- ggplot(list.sf[[i]])+geom_sf(aes(fill = .data[[x]], col = .data[[x]])) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = rev.pal[grep(x, diffs)] , aesthetics = c("fill","col"), name = paste(legends[grep(x, diffs)],"\n",names(list.sf)[i], sep = ""), limits = c(min.lim[[x]], max.lim[[x]])) +
# scale_fill_gradient2(low=palettes[3, grep(x,names(palettes))],mid=palettes[2, grep(x,names(palettes))],high = palettes[1, grep(x,names(palettes))],
# midpoint = 0, aesthetics = c("fill","col"),
# name = paste(legends[grep(x, diffs)],"\n",names(list.sf)[i], sep = ""), limits = c(min.lim[[x]], max.lim[[x]])) +
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = g, filename = paste("WR/", names(list.sf)[i], ".", x, ".png",sep = ""),dpi = "retina", width = 12, height = 16, units = "cm")
}
}
white <- ggplot()+theme_void()
ggsave(plot = white, filename = "white.png",dpi = "retina", width = 12, height = 16, units = "cm")listgrob <- list("WR/SSP126.2050.diffTemp.png", "WR/SSP126.2100.diffTemp.png","WR/SSP370.2050.diffTemp.png","WR/SSP370.2100.diffTemp.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
Figure 7: Map of forecasted temperature (ºC) for 2050 and 2100, under SSP 1-2.6 and SSP 3-7.0
listgrob <- list("WR/SSP126.2050.diffPrecip.png","WR/SSP126.2100.diffPrecip.png","WR/SSP370.2050.diffPrecip.png","WR/SSP370.2100.diffPrecip.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
Figure 12: Forecast of mean precipitation intensities for the period 2041-2060 (i.e., 2050) and for the period 2081-2100 (i.e., 2100), and change in precipitation compared to the average for 1995 to 2050, under SSP 1-2.6 (2ºC target) and SSP 3-7.0 (realistic worst-case) scenarios.
listgrob <- list("WR/SSP126.2050.diffNDVI.png","WR/SSP126.2100.diffNDVI.png","WR/SSP370.2050.diffNDVI.png","WR/SSP370.2100.diffNDVI.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
Figure 13: Forecast of mean NDVI for the period 2041-2060 (i.e., 2050) and for the period 2081-2100 (i.e., 2100), and change in NDVI compared to the average for 1995 to 2050, under SSP 1-2.6 (2ºC target) and SSP 3-7.0 (realistic worst-case) scenarios.
listgrob <- list("WR/SSP126.2050.diffRunoff.png","WR/SSP126.2100.diffRunoff.png","WR/SSP370.2050.diffRunoff.png","WR/SSP370.2100.diffRunoff.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
Figure 14: Forecast of mean runoff intensities for the period 2041-2060 (i.e., 2050) and for the period 2081-2100 (i.e., 2100), and change in runoff compared to the average for 1995 to 2050, under SSP 1-2.6 (2ºC target) and SSP 3-7.0 (realistic worst-case) scenarios.
listgrob <- list("WR/SSP126.2050.diffTNdep.png","WR/SSP126.2100.diffTNdep.png","WR/SSP370.2050.diffTNdep.png","WR/SSP370.2100.diffTNdep.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
Figure 15: Forecast of mean total nitrogen deposition for the period 2041-2060 (i.e., 2050) and for the period 2081-2100 (i.e., 2100), and change in TNdep compared to the average for 1995 to 2050, under SSP 1-2.6 (2ºC target) and SSP 3-7.0 (realistic worst-case) scenarios.
3.2 Forecast TOC for 2050 and 2100
The average concentration of total organic carbon in the coastal drainage basins was predicted based on the SELM with 5 predictors, fitted on the Fennoscandian dataset. The difference between the predicted log(TOC) for 2050 and 2100 under the SSP 1-2.6 and 3-7.0 scenarios, compared to the log(TOC) predicted for 1995, was computed and plotted.
fennoscandia <- readRDS("fennoscandia.33.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2050.rds")
temperature <- c(fennoscandia$Temp, sem.pred.wr$Temp,wr.ssp126.2050$Temp,wr.ssp370.2050$Temp, wr.ssp126.2100$Temp,wr.ssp370.2100$Temp)
precipitation <- c(fennoscandia$Precip, sem.pred.wr$Precip,wr.ssp126.2050$Precip,wr.ssp370.2050$Precip, wr.ssp126.2100$Precip,wr.ssp370.2100$Precip)
ndvi <- c(fennoscandia$NDVI, sem.pred.wr$NDVI,wr.ssp126.2050$NDVI,wr.ssp370.2050$NDVI, wr.ssp126.2100$NDVI,wr.ssp370.2100$NDVI)
runoff <- c(fennoscandia$Runoff, sem.pred.wr$Runoff,wr.ssp126.2050$Runoff,wr.ssp370.2050$Runoff, wr.ssp126.2100$Runoff,wr.ssp370.2100$Runoff)
tndep <- c(fennoscandia$TNdep, sem.pred.wr$TNdep,wr.ssp126.2050$TNdep,wr.ssp370.2050$TNdep, wr.ssp126.2100$TNdep,wr.ssp370.2100$TNdep)
toc <- c(fennoscandia$TOC, sem.pred.wr$TOC95,wr.ssp126.2050$TOC,wr.ssp370.2050$TOC, wr.ssp126.2100$TOC,wr.ssp370.2100$TOC)
neworder <- c("Fennoscandia","WR1995","SSP126.2050","SSP126.2100","SSP370.2050","SSP370.2100")
model <- factor(c(rep("Fennoscandia",dim(fennoscandia)[1]),
rep("WR1995",dim(sem.pred.wr)[1]),
rep("SSP126.2050", dim(wr.ssp126.2050)[1]),
rep("SSP370.2050", dim(wr.ssp370.2050)[1]),
rep("SSP126.2100", dim(wr.ssp126.2100)[1]),
rep("SSP370.2100", dim(wr.ssp370.2100)[1])),
levels = neworder)
compare.ssp <- data.frame(value = c(temperature, precipitation, runoff, ndvi, toc), variable = c(rep("Temperature (°C)", length(temperature)), rep("Precipitation (mm)", length(precipitation)), rep("Runoff (mm/y)",length(runoff)), rep("NDVI",length(ndvi)), rep("TOC (mg/L)", length(toc))), model = rep(model,5))
ggplot(compare.ssp)+geom_boxplot(aes(y= value,col=variable))+
labs(y="", x = "")+scale_color_viridis_d(end = 0.8)+
theme_minimal()+
theme(legend.position = "", axis.text.x = element_blank())+
facet_grid(cols = vars(model), rows = vars(factor(variable, levels = c("Temperature (°C)","Precipitation (mm)","Runoff (mm/y)","NDVI","TOC (mg/L)"))), scales = "free_y", shrink = T)
ggsave("WR/compare.ssp.png", width = 30, height = 25, units = "cm", dpi = "retina")knitr::include_graphics("WR/compare.ssp.png")
Figure 16: Boxplots of the main predictors for the fennoscandian dataset, the coastal drainage basins in 1995, and the forecasts for the coastal drainage basins for SSP 1-2.6 and 3-7.0 in 2050 and 2100
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp126 <- predict(sem.fennoscandia.5,wr.sf.126.2050,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp126 <- cbind(wr.sf.126.2050,wr.pred.ssp126)
wr.ssp126$TOC <- exp(wr.ssp126$fit)
wr.ssp126 <- merge(wr.ssp126, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp126,"WR/wr.ssp126.2050.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp370 <- predict(sem.fennoscandia.5,wr.sf.370.2050,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp370 <- cbind(wr.sf.370.2050,wr.pred.ssp370)
wr.ssp370$TOC <- exp(wr.ssp370$fit)
wr.ssp370 <- merge(wr.ssp370, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp370,"WR/wr.ssp370.2050.rds")sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp126 <- predict(sem.fennoscandia.5,wr.sf.126.2100,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp126 <- cbind(wr.sf.126.2100,wr.pred.ssp126)
wr.ssp126$TOC <- exp(wr.ssp126$fit)
wr.ssp126 <- merge(wr.ssp126, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp126,"WR/wr.ssp126.2100.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp370 <- predict(sem.fennoscandia.5,wr.sf.370.2100,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp370 <- cbind(wr.sf.370.2100,wr.pred.ssp370)
wr.ssp370$TOC <- exp(wr.ssp370$fit)
wr.ssp370 <- merge(wr.ssp370, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp370,"WR/wr.ssp370.2100.rds")The average TOC concentration was plotted.
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2100.rds")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2100.rds")
list.sf <- list(wr.ssp126.2050, wr.ssp126.2100, wr.ssp370.2050, wr.ssp370.2100)
names(list.sf) <- c("SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
for(i in names(list.sf)){
list.sf[[i]]$diffTOC <- list.sf[[i]]$TOC - list.sf[[i]]$TOC95
list.sf[[i]]$diff.log <- list.sf[[i]]$fit - list.sf[[i]]$wr.pred.fit
list.sf[[i]]$exp.diff.log <- exp(list.sf[[i]]$diff.log)
list.sf[[i]]$percent.change <- (list.sf[[i]]$exp.diff.log - 1) * 100
}
all.diff <- data.frame(gid = c(wr.ssp126.2050$gid, wr.ssp370.2050$gid, wr.ssp126.2100$gid, wr.ssp370.2100$gid))
for(x in c("fit","diffTOC","diff.log","exp.diff.log","percent.change")){
all.vec <- c()
for(i in names(list.sf)){
vec <- dplyr::select(st_drop_geometry(list.sf[[i]]),all_of(x))
all.vec <<- rbind(all.vec,vec)
}
all.diff <- cbind(all.diff, all.vec)
}
max.lim <- apply(all.diff, 2, max, na.rm = T)
min.lim <- apply(all.diff, 2, min, na.rm = T)
t0 <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=wr.pred.fit, col= wr.pred.fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "logTOC in 1995")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/toc_1995.png", t0 ,dpi = "retina", width = 12, height = 16, units = "cm")
t1 <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "a) Forecast logTOC\n SSP1-2.6 - 2050")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
# scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Forecast logTOC\n SSP1-2.6 - 2050", mid = (max.lim[2]-min.lim[2])/2, limits = c(min.lim[2],max.lim[2]), rev = T) +
# theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp126.2050.png", t1 ,dpi = "retina", width = 12, height = 16, units = "cm")
t2 <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "b) Forecast logTOC\n SSP1-2.6 - 2100")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp126.2100.png", t2 ,dpi = "retina", width = 12, height = 16, units = "cm")
t3 <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "c) Forecast logTOC\n SSP3-7.0 - 2050")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp370.2050.png", t3 ,dpi = "retina", width = 12, height = 16, units = "cm")
t4 <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "d) Forecast logTOC\n SSP3-7.0 - 2100")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp370.2100.png", t4 ,dpi = "retina", width = 12, height = 16, units = "cm")
knitr::include_graphics("WR/toc_1995.png")
list.grob <- c("WR/logTOC_ssp126.2050.png", "WR/logTOC_ssp126.2100.png","WR/logTOC_ssp370.2050.png","WR/logTOC_ssp370.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow = 2, ncol = 2, labels = c("a)","b)","c)","d)"), label_size = 10)
Figure 17: Forecast of mean log(TOC) for the period 2041-2060 (i.e., 2050) and for the period 2081-2100 (i.e., 2100), under SSP 1-2.6 (2ºC target) and SSP 3-7.0 (realistic worst-case) scenarios.
p1 <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "a) Forecasted percent increase in [TOC]\n SSP1-2.6 - 2050")+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Forecast [TOC] increase in %\n SSP1-2.6 - 2050", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
p2 <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="b) Forecast [TOC] increase in %\n SSP1-2.6 - 2100", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "b) Forecasted percent change TOC\n SSP1-2.6 - 2100")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
p3 <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="c) Forecast [TOC] increase in %\n SSP3-7.0 - 2050", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "c) Forecasted percent change TOC\n SSP3-7.0 - 2050")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
p4 <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="d) Forecast [TOC] increase in %\n SSP3-7.0 - 2100", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "d) Forecasted percent change TOC\n SSP3-7.0 - 2100")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
pp <- grid.arrange(p1, p2, p3, p4, ncol = 2, nrow = 2)
ggsave("forecast.percent.change.toc.png", pp, width = 12, height = 16, units = "cm")
highchangeid <- subset(list.sf$SSP370.2100, exp.diff.log > 25)$gid
pp1 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2050"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "a) Forecasted percent change TOC\n SSP1-2.6 - 2050")+
geom_sf(data = subset(list.sf[["SSP126.2050"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
pp2 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2100"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "b) Forecasted percent change TOC\n SSP1-2.6 - 2100")+
geom_sf(data = subset(list.sf[["SSP126.2100"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
pp3 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2050"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "c) Forecasted percent change TOC\n SSP3-7.0 - 2050")+
geom_sf(data = subset(list.sf[["SSP370.2050"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
pp4 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2100"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "d) Forecasted percent change TOC\n SSP3-7.0 - 2100")+
geom_sf(data = subset(list.sf[["SSP370.2100"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ppp <- grid.arrange(pp1, pp2, pp3, pp4, ncol = 2, nrow = 2)
ggsave("WR/forecast.prop.change.png", ppp, width = 12, height = 16)
c1 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2050"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Forecast [TOC] change in %\n SSP1-2.6 - 2050", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP126.2050"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp126.2050.png", c1 ,dpi = "retina", width = 12, height = 16, units = "cm")
c2 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2100"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="b) Forecast [TOC] change in %\n SSP1-2.6 - 2100", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP126.2100"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp126.2100.png", c2 ,dpi = "retina", width = 12, height = 16, units = "cm")
c3 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2050"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="c) Forecast [TOC] change in %\n SSP3-7.0 - 2050", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP370.2050"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp370.2050.png", c3 ,dpi = "retina", width = 12, height = 16, units = "cm")
c4 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2100"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="d) Forecast [TOC] change in %\n SSP3-7.0 - 2100", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP370.2100"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp370.2100.png", c4 ,dpi = "retina", width = 12, height = 16, units = "cm")
list.grob <- c("WR/percentTOC_ssp126.2050.png", "WR/percentTOC_ssp126.2100.png","WR/percentTOC_ssp370.2050.png","WR/percentTOC_ssp370.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow = 2, ncol = 2, labels = c("a)","b)","c)","d)"), label_size = 10)
Figure 18: Forecast of change of TOC concentration in % for the period 2041-2060 (i.e., 2050) and for the period 2081-2100 (i.e., 2100), under SSP 1-2.6 (2ºC target) and SSP 3-7.0 (realistic worst-case) scenarios. In black: catchments with concentration change superior to 307% (5% of the catchments)
g <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Difference logTOC\n SSP1-2.6 - 2050", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2050.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="b) Difference logTOC\n SSP1-2.6 - 2100", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2100.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="c) Difference logTOC\n SSP3-7.0 - 2050", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2050.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="d) Difference logTOC\n SSP3-7.0 - 2100", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2100.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")listgrob <- list("WR/SSP126.2050.diff.log.png", "WR/SSP126.2100.diff.log.png","WR/SSP370.2050.diff.log.png","WR/SSP370.2100.diff.log.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
ssp <- plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
save_plot("ssp.toc.png",ssp, ncol = 2, nrow = 2, base_height = 4, base_width = 3)
print(ssp)
for(i in names(list.sf)){
list.sf[[i]]$sign.diff <- sign(list.sf[[i]]$diffTOC) %>% as.factor()
}
g <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=sign.diff, col= sign.diff), size= 0.05)+
scale_color_manual(values = c("lightblue","firebrick"), labels = c("Decreased [TOC]", "Increased [TOC]"), aesthetics = c("fill","col"), name = "SSP 1-2.6 2050")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2050.sign.diff.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=sign.diff, col= sign.diff), size= 0.05)+
scale_color_manual(values = c("lightblue","firebrick"), labels = c("Decreased [TOC]", "Increased [TOC]"), aesthetics = c("fill","col"), name = "SSP 1-2.6 2100")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2100.sign.diff.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=sign.diff, col= sign.diff), size= 0.05)+
scale_color_manual(values = c("lightblue","firebrick"), labels = c("Decreased [TOC]", "Increased [TOC]"), aesthetics = c("fill","col"), name = "SSP 3-7.0 2050")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2050.sign.diff.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=sign.diff, col= sign.diff), size= 0.05)+
scale_color_manual(values = c("lightblue","firebrick"), labels = c("Decreased [TOC]", "Increased [TOC]"), aesthetics = c("fill","col"), name = "SSP 3-7.0 2100")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2100.sign.diff.png",dpi = "retina", width = 12, height = 16, units = "cm")listgrob <- list("WR/SSP126.2050.sign.diff.png", "WR/SSP126.2100.sign.diff.png","WR/SSP370.2050.sign.diff.png","WR/SSP370.2100.sign.diff.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
sdf <- plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
save_plot("sign.diff.toc.png",sdf, ncol = 2, nrow = 2, base_height = 4, base_width = 3)
print(sdf)
4 Coastal drainage basins: coastal darkening
Most coastal drainage basins in Sweden and Finland drain into to the Baltic Sea (see Supplementary 2), while Norwegian sea watershed drain in Skagerrak, in the North sea, in the Norwegian Sea and in the Barent sea in the North (Sætre 2007). The drainage basin of the Baltic sea matches almost exactly the borders of Sweden and Finland : https://www.grida.no/resources/5324. All contribute to the Norwegian Coastal Current water masses. The estimated exported TOC was computed as:
\[ TOC_exp = [TOC] (mg.L^{-1}) \times Runoff (mm/y, or L.m^{-2}.y^{-1}) \times Watershed area (m^2)\]
Where the TOC concentration (past and future) was obtained by back-transforming the results of the model (log(TOC)) into TOC in mg/L. The result of the above equation, giving the amount of TOC exported in mg, is converted to Tg. Forecast are presented on Figure 19*.
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.centroid <- st_centroid(wr.sf.95, of_largest_polygon = T) %>% dplyr::select("gid","geometry")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
nor <- st_read("Country_shapefile/norway.shp") %>% st_transform(st_crs(wr.centroid))
swe <- st_read("Country_shapefile/sweden.shp") %>% st_transform(st_crs(wr.centroid))
fin <- st_read("Country_shapefile/finland.shp") %>% st_transform(st_crs(wr.centroid))
wr.norge <- wr.centroid[nor,]
wr.norge$Nation <- "Norway"
wr.sweden <- wr.centroid[swe,]
wr.sweden$Nation <- "Sweden"
wr.finland <- wr.centroid[fin,]
wr.finland$Nation <- "Finland"
wr.nation <- rbind(wr.norge,wr.sweden,wr.finland)
wr.sf.95 <- wr.sf.95 %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid") %>% sp::merge(dplyr::select(st_drop_geometry(sem.pred.wr), c("gid","wr.pred.fit","TOC95")), by = "gid")
names(wr.sf.95[which(names(wr.sf.95) == "TOC95")]) <- "TOC"
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")%>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2100.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2100.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
list.sf <- list(wr.sf.95, wr.ssp126.2050, wr.ssp126.2100, wr.ssp370.2050, wr.ssp370.2100)
names(list.sf) <- c("Data1995","SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
export <- function(sf){
discharge <- sf$Runoff*sf$area # Q in L/y
toc.export <- discharge * sf$TOC * 10^(-15) # amount of exported TOC in Tg/y
sf2 <- cbind(sf,discharge,toc.export)
return(sf2)
}
list.sf <- lapply(list.sf, export)
list.mean <- c()
for(i in 1:length(list.sf)){
list.mean[[i]] <- list.sf[[i]] %>% st_drop_geometry() %>% group_by(Nation) %>% summarise(total = sum(toc.export)) %>% as.data.frame()
}
names(list.mean) <- names(list.sf)
list.diff <- c()
for(i in 1:length(list.mean)){
df1 <- list.mean[[1]]
df2 <- list.mean[[i]]
diff <- df2[,2] - df1[,2]
x <- cbind(list.mean[[i]],diff)
names(x) <- c("Nation", names(list.mean)[i], paste("diff", names(list.mean)[i], sep = ""))
list.diff[[i]] <- x
}
toc.export.df <- list.diff[[1]] %>% merge(list.diff[[2]], by = "Nation") %>% merge(list.diff[[3]], by = "Nation") %>% merge(list.diff[[4]], by = "Nation") %>% merge(list.diff[[5]], by = "Nation") %>% melt(id.vars = c("Nation"), value.name = "TOC_Tg")
toc.export.df$Model <- NA
m.126 <- grep("126",toc.export.df$variable)
m.370 <- grep("370",toc.export.df$variable)
toc.export.df$Model[m.126] <- "SSP126"
toc.export.df$Model[m.370] <- "SSP370"
toc.export.df$Model[-c(m.126,m.370)] <- "Historical"
toc.export.df$Year <- NA
m.95 <- grep("95",toc.export.df$variable)
m.2050 <- grep("2050",toc.export.df$variable)
m.2100 <- grep("2100",toc.export.df$variable)
toc.export.df$Year[m.95] <- "1995"
toc.export.df$Year[m.2050] <- "2050"
toc.export.df$Year[m.2100] <- "2100"
toc.export.df$Diff <- NA
m.d <- grep("diff",toc.export.df$variable)
toc.export.df$Diff[m.d] <- "Difference"
toc.export.df$Diff[-m.d] <- "Absolute"
toc.export.sum <- toc.export.df %>% group_by(Model,Year,Diff) %>% summarise(TOC_Tg = sum(TOC_Tg))
toc.export.sum$Nation <- "Total"
toc.export.df <- rbind(toc.export.sum,toc.export.df)
toc.export.df <- transform(toc.export.df, Nation = factor(Nation, levels = c("Norway","Sweden","Finland","Total")))
toc.export.df$variable <- NULL
toc.export.df$TOC_Tg <- as.numeric(toc.export.df$TOC_Tg)
knitr::kable(toc.export.df) %>% kable_styling()
saveRDS(toc.export.df,"WR/toc.export.df.rds")
ggplot() +
geom_point(data = toc.export.df, aes(x=Year, y = TOC_Tg, col = Model, shape = Model)) +
geom_path(data = toc.export.df, aes(x=Year, y = TOC_Tg, col = Model, group = Model), size = 0.5) +
geom_abline(slope = 0, intercept = 0, col = "firebrick4", size = 0.3)+
scale_color_manual(values = c("firebrick4","steelblue4","darkorange"))+
geom_label(data = filter(toc.export.df, Model == "SSP126" & Diff == "Absolute"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = -0.5, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "SSP126"& Diff == "Difference"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = -0.1, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "SSP370"& Diff == "Absolute"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.5, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "SSP370"& Diff == "Difference"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.1, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "Historical"& Diff == "Absolute"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.5, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "Historical"& Diff == "Difference"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.1, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
labs(x="Year", y= "TOC (Tg)")+
facet_grid(cols=vars(Nation), rows = vars(Diff), scales = "free_y")+
theme_bw(base_size = 12)+
theme(legend.position = "bottom")
ggsave(filename = "WR/TOC_export.png", dpi = "retina", width = 20, height = 15, units = "cm")
Table 1: TOC exported from coastal drainage basins (Tg)
toc.export.df <- readRDS("WR/toc.export.df.rds") %>% dplyr::select(c("Nation", "Model", "Year","Diff","TOC_Tg"))
toc.export.df %>% arrange(Model) %>%
knitr::kable(digits = 3,
caption = "TOC export by country",
col.names = c("Country","Model","Year","Absolute/differential export","TOC export (Tg)")) %>%
kable_styling(bootstrap_options = "bordered",
position = "center",
full_width = F)
knitr::include_graphics("WR/TOC_export.png")
Figure 20: TOC exported from coastal drainage basins (Tg)
References
CORDEX. 2021. “CORDEX experiment design for dynamical downscaling of CMIP6,” no. May: 1–8. https://cordex.org/wp-content/uploads/2021/05/CORDEX-CMIP6_exp_design_RCM.pdf.
Detsch, Florian. 2021. “Download and Process GIMMS NDVI3g Data [R package gimms version 1.2.1].” Comprehensive R Archive Network (CRAN). https://cran.r-project.org/package=gimms.
European Environment Agency. 2021. “Emissions of the main air pollutants in Europe.” https://www.eea.europa.eu/ims/emissions-of-the-main-air.
Finstad, Anders G. 2017. “Environmental data NINAnor/NOFA Wiki.” https://github.com/NINAnor/NOFA/wiki/Environmental-data.
Hausfather, Zeke. 2019. “CMIP6: the next generation of climate models explained.” https://www.carbonbrief.org/cmip6-the-next-generation-of-climate-models-explained.
Hijmans, Robert J, Jacob Van Etten, Joe Cheng, Michael Sumner, Matteo Mattiuzzi, Jonathan A Greenberg, Oscar Perpinan, et al. 2018. “Package ’raster’ Type Package Title Geographic Data Analysis and Modeling.” https://github.com/rspatial/raster/issues/ https://cran.r-project.org/package=raster.
Masson-Delmotte, V., P. Zhai, A. Pirani, S. L. Connors, C. Péan, S. Berger, N. Caud, et al. 2021. “IPCC. 2021: Summary for Policymakers.” In Climate Change 2021: The Physical Science Basis. Contribution of Working Group i to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge, 3–32. Cambridge, United Kingdom; New Yoak, USA: IPCC. https://www.ipcc.ch/report/ar6/wg1/.
Myers-Smith, Isla H, Jeffrey T Kerby, Gareth K Phoenix, Jarle W Bjerke, Howard E Epstein, Jakob J Assmann, Christian John, et al. 2020. “Complexity revealed in the greening of the Arctic.” https://doi.org/10.1038/s41558-019-0688-1.
Raju, Komaragiri Srinivasa, and Dasika Nagesh Kumar. 2020. “Review of approaches for selection and ensembling of GCMS.” Journal of Water and Climate Change 11 (3): 577–99. https://doi.org/10.2166/wcc.2020.128.
Sætre, Roald. 2007. The Norwegian Coastal Current - Oceanography and Climate. Edited by Tapir Academic Publisher. Trondheim, Norway.
The National Center for Atmospheric Research. 2018. “Global GIMMS NDVI3g v1 dataset (1981-2015).” A Big Earth Data Platform for Three Poles. http://poles.tpdc.ac.cn/en/data/9775f2b4-7370-4e5e-a537-3482c9a83d88/.
Tucker, Compton J, Jorge E Pinzon, Molly E Brown, Daniel A Slayback, Edwin W Pak, Robert Mahoney, Eric F Vermote, and Nazmi El Saleous. 2005. “An extended AVHRR 8-km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data.” International Journal of Remote Sensing 26 (20): 4485–98. https://doi.org/10.1080/01431160500168686.